/*


@author Vishwajeet Dusan
*/


#include"matrix.h"

class LUDecomposition:public Matrix
{
  public :
    
      Matrix Li,Ui,A,B,Xn,Yn;   
   
      
   // calculate Li.matrix= aij - (Li.matrixkUkj ) where k<(k-1) and i<j        
   void calculateLi(int,int);
   
   // calculate Ui.matrix= (aij - (Li.matrixkUkj ))/Li.matrixi where k<(k-1) and i>=j        
   void calculateUi(int ,int);
   
   //intitial value of Li.matrix
   void initialLi();
   
   //initial value of Ui.matrix
   void initialUi();
   
   //test whether the calculated Ui.matrix and Li.matrix are correct
   //  By Li.matrix*Ui.matrix = A
   void testLUD();
   
    
   //Collection  of calculateLi ,CalculateUi and other method's 
   //pass parameter input example file which contains matrix A and matrix B
    void luDeco(char [50],char [50]);
   
   Matrix findRoots();
   
};

 Matrix LUDecomposition::findRoots()
 {
 
  
  initialLi();
  
  initialUi();
  
  for(int i=1;i<A.row;i++)
   {
    for(int j=1;j<A.column;j++)
       {
       
             if(i>=j)
	     {   
	       	
	//	cout<<"out Li i= "<<i<<"  j= "<<j<<endl;
		calculateLi(i,j);
	     }	
	     else	
	      {
	        
	//	cout<<"out Ui i= "<<i<<"  j= "<<j<<endl;
		calculateUi(i,j);    
	      }	
		
       }
   
   }   
 
  //Li.printMatrix();
 // Ui.printMatrix();

 // cout<<"After Printing LiUi"<<endl<<endl;
  
  
  Yn=Li.getValueOfEquation(B);
    
   //Yn.printMatrix();
  
  //cout<<"After Printing Yn"<<endl<<endl; 
  Xn=Ui.getValueOfEquation(Yn);
  
   cout<<"Solution for augmented matrix is\n";
   Xn.printMatrix();
  

return Xn; 
 }
 
 
void LUDecomposition::luDeco(char example[50],char exampleB[50])
{
  
   A.readMatrix(example);
   B.readMatrix(exampleB);
   
   Li.setMatrixSize(A.row,A.column);
   Li.makeElement(0);
      
   Ui.setMatrixSize(A.row,A.column);
   Ui.makeElement(0);
  
    findRoots(); 
  
 

}


void LUDecomposition::initialUi()
{
   
   for(int i=0;i<A.column;i++)
   {
     Ui.matrix[0][i]=A.matrix[0][i]/(double)A.matrix[0][0]; 
     Ui.matrix[i][i]=1;
   }
   
}

void LUDecomposition::initialLi()
{

 for(int i=0;i<A.row;i++)
 {
    Li.matrix[i][0]=A.matrix[i][0]; 
 }
 

}

void LUDecomposition::calculateUi(int i,int j)
{
   double U;
        U=0;
	
	
	
        for(int k=0;k<j;k++)
	    {
	        U=U+(Li.matrix[i][k]*Ui.matrix[k][j]);    
	    
	    }
	   Ui.matrix[i][j]=A.matrix[i][j]-U; 
	   Ui.matrix[i][j]=Ui.matrix[i][j]/Li.matrix[i][i]; 
	
   
}


void LUDecomposition::calculateLi(int i,int j)
{
double L;
  L=0;
   	   for(int k=0;k<j;k++)
	    {
	     L=L+(Li.matrix[i][k]*Ui.matrix[k][j]);
	     
	  //   cout<<"In Li L= "<<L;
	   //  cout<<"In Li.matrix[i][k]*Ui.matrix[k][j]= "<<Li.matrix[i][k]*Ui.matrix[k][j]<<endl;
	    }
	   Li.matrix[i][j]=A.matrix[i][j]-L;
	   //cout<<"Li  "<<Li.matrix[i][j]<<" A"<<A.matrix[i][j]<<endl<<endl;
	
   
}


void LUDecomposition::testLUD()
{
  
   Matrix temp(A.row,A.column);
   
     temp=Li*Ui;
     if(temp==A)
      cout<<"Li*Ui==A";
     else
     {
      cout<<"Li*Ui!=A";  
      //temp.printMatrix();
      }
  
}
